\documentclass[12pt]{article} \usepackage{amsmath,amsfonts}
%\usepackage{showlabels}
\usepackage[pdftex]{graphicx,color}

\newcommand{\htop}{{h_{\text{top}}}}
\newcommand{\T}{{\cal P}}
\newcommand{\logprob}{{\cal L}}
\newcommand{\measure}{\mu}
\newcommand{\Aop}{{\cal A}}
\newcommand{\Aindicate}{\alpha}
\newcommand{\slope}{s}
\newcommand{\density}{p}
\newcommand{\stationary}{\mu}
\newcommand{\function}{\phi}
\newcommand{\rightfunction}{{_{_{\scalebox{.4}{R}}} \rho}}
\newcommand{\leftfunction}{{_{_{\scalebox{.4}{L}}} \rho}}
\renewcommand{\th}{^{\text{th}}}
\newcommand{\field}[1]{\mathbb{#1}}
\newcommand\REAL{\field{R}}

\title{Uniform Probability for a Continuous Time Monotonic Process }

\author{Andrew M.\ Fraser}
\begin{document}
\maketitle
\begin{abstract}
  I describe a probability measure for functions of time that are
  constrained both to lie in within certain bounds and to be
  monotonic.  In coordinates that make the constraints invariant with
  respect to shifts in time, the measure corresponds to a stationary
  stochastic process.  Roughly, the measure makes all allowed
  sequences have the same probability with respect to a norm.
\end{abstract}

\section*{To Do:}
\label{sec:do}

\begin{itemize}
\item Understand behavior of eigenvalues and eigenfunctions in
  $\tau \rightarrow 0$ limit.  Include results of numerical
  calculations (with $n_g=3000$) in plots of eigenvalues and
  eigenvectors.
\item Figure~\ref{fig:paths}: Plot of
  $\tilde T_{t=\frac{1}{4},g_0=\frac{1}{2},n}$ for several values of
  $n$.
\item Surface plot of $p(g(t_1),g(t_0))$ for $t_0=0$ and
  $t_t=\frac{1}{4}$.
\end{itemize}

\section{Introduction}
\label{sec:introduction}
\begin{description}
\item[Motivation:] Priors for regions in function space constrained by
  physical laws.  Cite Hixson2000, Fritz1996,
  Pemberton2011LA-UR-11-04999, and FraserLA-UR-13-21824
\item[Connection to Information Theory:] Theorem 8 of Shannon.  Define
  stationary and Markov for discrete time discrete value process.
\item[Connection to Statistical Mechanics:] Jaynes maximum entropy
  derivation of Maxwell Boltzmann distribution.  Coordinate dependent.
\item[Connection to path integrals:]
\item[Not Levy process:] Connection to Brownian motion
\item[Outline:]
\end{description}

\section{Specific Problem}
\label{sec:specific}

I seek a \emph{uniform} distribution on a \emph{restricted} set of
functions.  Figure~\ref{fig:bounds} illustrate the set which I define
by the following constraints:
\begin{subequations}
  \label{eq:f_def}
  \begin{description}
  \item[Bounded:] For all values of the independent variable $t$ the
    function lies between $1-t$ and $-t$,
    \begin{equation}
      \label{eq:bounded_f}
      -t \leq f(t) \leq 1-t
    \end{equation}
  \item[Monotonic:] The function is monotonic,
    \begin{equation}
      \label{eq:monotonic_f}
      f(t+\delta) \leq f(t) ~ \forall (t,\delta): \delta \geq 0
    \end{equation}
  \end{description}
\end{subequations}

\begin{figure*}
  \centering
  \resizebox{0.75\textwidth}{!}{\includegraphics{bounds.pdf}}
  \caption{Bounds on allowed functions}
  \label{fig:bounds}
\end{figure*}
If one derives a new function $g$ by shifting $f$ by $t$, ie,
\begin{equation}
  \label{eq:shift}
  g(t) = f(t) + t,
\end{equation}
then the constraints \eqref{eq:f_def} on $f$ imply the following
constraints on $g$:
\begin{subequations}
  \label{eq:g_def}
  \begin{align}
    \label{eq:bounded_g}
    0 &\leq g(t) \leq 1 ~ \forall t\\
    \label{eq:monotonic_g}
    g(t+\delta) &\leq g(t)+\delta ~ \forall(t,\delta): \delta \geq 0  
  \end{align}
\end{subequations}

While Eqns.~\eqref{eq:g_def} describe the constraints that pretty
simply, the goal of this paper is the more difficult task of defining
a uniform distribution over the set.  I require that the probability
measure $\measure$ have the following properties:
\begin{description}
\item[Stationary:] 
\item[Consistent:] 
\end{description}


\section{Continuous $g$:}
\label{sec:continuous_g}

I let $\Aop_\tau$ denote the operator for $\tau = \slope \Delta t$
with continuous $g$.  (Note: $\tau$ is the only parameter in the
problem.)  In turn I define $\Aop_\tau$ in terms of the function
$\Aindicate_\tau$ that has the value one if its arguments are an
allowed pair $\left( g(0), g(\tau) \right)$ and zero otherwise, viz:
\begin{subequations}
  \begin{align}
    \label{eq:indicate}
    \Aindicate(u,v) &\equiv
    \begin{cases}
      1 & 0 \leq u \leq 1,~ 0 \leq v \leq 1,~ v \leq u + \tau\\
      0 & \text{otherwise}
    \end{cases}\\
    \label{eq:operateright}
    \left( \Aop_\tau  \function \right)(g) &\equiv \int
            \Aindicate(g,h) \function(h) dh \\
    \label{eq:operateleft}
    \left( \function \Aop_\tau \right)(g) &\equiv \int
           \function(h) \Aindicate(h,g) dh \\
    &= \int_{\max(g-\tau,0)}^1 \function(h) dh.
  \end{align}
\end{subequations}
I use the following notation for an eigenvalue of $\Aop_\tau$ and
corresponding \emph{left} and \emph{right} eigenfunctions:
\begin{align*}
 \left( \leftfunction_\lambda \Aop_\tau \right)(g) &= \lambda
 \leftfunction_\lambda (g)\\
 \left(\Aop_\tau  \rightfunction_\lambda \right)(g) &= \lambda
 \rightfunction_\lambda (g).
\end{align*}
A symmetry argument yields
\begin{equation}
  \label{eq:symmetry}
  \rightfunction_\lambda(g) \propto \leftfunction_\lambda(1-g).
\end{equation}
In the following, since I am only interested in the largest
eigenvalue and the corresponding eigenfunctions\footnote{Called
  the Peron Frobenius eigenvalue and eigenvectors}, I either don't use
subscripts or use them to denote $\tau=s\Delta t$.

Starting from this definition of an eigenfunction
\begin{equation}
  \label{eq:1}
  \lambda_\tau \leftfunction_\tau(g) = \left( \leftfunction_\tau
    \Aop_\tau \right) (g) = \int_0^1 \leftfunction_\tau(h)
  \Aindicate_\tau(h,g) dh,  
\end{equation}
I derive a procedure for constructing eigenfunctions.  Note that
$\Aindicate_\tau(h,g) = 1$ for all $h$ if $g \leq \tau$.
Consequently, the integral in \eqref{eq:1} has the same value for any
$g$ less than $\tau$, ie,
\begin{equation*}
  \int_0^1 \leftfunction_\tau(h) dh \equiv a ~\forall g: g \leq \tau.
\end{equation*}
I can choose to scale the eigenfunction so that $a=\lambda_\tau$ and
\begin{equation*}
  \leftfunction_\tau  (g) = 1 ~\forall g: g \leq \tau.
\end{equation*}
For larger values of $g$ the integrand in \eqref{eq:1} is one for
$g-\tau \leq h < 1$ and
\begin{align}
  \lambda_\tau \leftfunction_\tau(g) &= \int_{g-\tau}^1
      \leftfunction_\tau(h) dh ~~\forall g > \tau \nonumber \\
  &= \int_0^1 \leftfunction_\tau(h) dh - \int_0^{g-\tau}
    \leftfunction_\tau(h) dh \nonumber \\
  &= \lambda_\tau - \int_0^{g-\tau} \leftfunction_\tau(h) dh \nonumber
  \\
  \label{eq:delay_ode}
  \frac{\partial}{\partial g} \leftfunction_\tau(g) &=
            \frac{-1}{\lambda_\tau} \leftfunction_\tau(g-\tau)
             ~~\forall g > \tau .  
\end{align}
Note that for $\tau = 0$, I get the following pair of equations that
has no solutions:
\begin{align*}
  \frac{\partial}{\partial g} \leftfunction_0(g) &=
       \frac{-1}{\lambda_0} \leftfunction_0(g) \\
  \lambda_0 &= \int_0^1 \leftfunction_0(h) dh.
\end{align*}

If $\tau = 1$ then the image of every function $\function$ is a
multiple of the constant function $\leftfunction_\tau(g) = 1$ and the
largest eigenvalue is $1$.  If $\tau = \frac{1}{2}$, then
\begin{equation*}
  \leftfunction_\tau(g) =
  \begin{cases}
    1 & g \leq \frac{1}{2} \\
    1 + \frac{1}{\lambda}\left( \frac{1}{2} - g \right) & g > \frac{1}{2}
  \end{cases},
\end{equation*}
and
\begin{align*}
  \lambda_\tau \leftfunction_\tau (g=0) &= \left( \leftfunction_\tau
                                          \Aop_\tau \right) (g=0) \\ 
  \lambda_\tau &= \left( \int_0^{\frac{1}{2}} dh + \int_{\frac{1}{2}}^1
                1 + \frac{1}{\lambda_\tau}\left( \frac{1}{2} - h \right)
                dh \right) \\
  &= 1 - \frac{1}{8 \lambda_\tau} \\
              &= \frac{1}{2} + \frac{1}{\sqrt{8}}.
\end{align*}

\subsection{Recursion Formula}
\label{sec:recursion}

If $\tau = \frac{1}{n}$, then $\leftfunction_\tau$ is piecewise
polynomial, and using $k$ as the index of the segments starting with
$k=0$, the highest power of $g$ in the $k\th$ segment is $k$.  By
defining $p_k$, a version of the $k\th$ segment rescaled so that it
starts at one, I derive a simple recursion formula for
$\leftfunction_\tau$.  With the definition
\begin{equation}
  \label{eq:pk_def}
  p_k(x) \equiv
  \begin{cases}
    \frac{\leftfunction_\tau(x+k*\tau)} {\leftfunction_\tau(k*\tau)} &
    0 \leq x \leq \tau \\
    \text{undefined} & \text{otherwise}
  \end{cases},
\end{equation}
I find
\begin{align*}
  p_0(x) &= 1 \\
  p_{k+1}(x) &= 1 - \frac{1}{\lambda} \int_0^{x} p_k(s) ds \\
  p_1(x) &= 1 - \frac{x}{\lambda} \\
  p_2(x) &= 1 - \frac{x}{\lambda} + \frac{x^2}{2\lambda^2} \\
  p_k(x) &= \sum_{n=0}^k \frac{1}{n!} \left( -\frac{x}{\lambda}
           \right)^n,
\end{align*}
and in the interval $k\tau \leq g \leq (k+1)\tau$ the eigenfunction is
\begin{equation}
  \label{eq:eigenfunction}
  \leftfunction_\tau(g) = p_k(g-k\tau) \prod_{n=0}^k p_n(\tau).
\end{equation}

\subsection{Approximate Limits}
\label{sec:limits}

Here I seek the behavior of $\lambda_\tau$ and $\leftfunction_\tau$ as
$\tau \rightarrow 0$.  Based on \eqref{eq:delay_ode}, I assume that
for small $\tau$ the eigenfunction has the following approximate form
\begin{equation}
  \label{eq:limit_function}
  \leftfunction(g) \approx
  \begin{cases}
    1 & g \leq \tau \\
    e^{- \frac{g - \tau}{\lambda}} & g > \tau .
  \end{cases}
\end{equation}
From \eqref{eq:1} in the region $g<\tau$ I write
\begin{align}
  \lambda &\approx \tau + e^{\frac{\tau}{\lambda}} \int_\tau^1
  e^{\frac{-h}{\lambda}} dh \nonumber  \\
  &= \tau - \lambda e^{\frac{\tau}{\lambda}} \left[
    e^{\frac{-h}{\lambda}} \right]_\tau^1\nonumber  \\
  &= \tau + \lambda \left( 1 - e^{\frac{\tau -1 }{\lambda}} \right)
    \nonumber \\
  \label{eq:limit_value}
  \tau & \approx \lambda e^{\frac{\tau - 1}{\lambda}}.
\end{align}
I let $\tilde \lambda(\tau)$ denote the value of $\lambda$ that solves
\eqref{eq:limit_value} as an equality.

Figures~\ref{fig:eigenfunctions} and~\ref{fig:eigenvalues} illustrate
$\leftfunction_{\frac{1}{n}}$ and $\lambda_{\frac{1}{n}}$ respectively
for a few values of $n$.\marginpar{Explain algorithms}

\begin{figure*}
  \centering
  \resizebox{0.75\textwidth}{!}{\includegraphics{eigenfunctions.pdf}}
  \caption{Eigenfunctions of $\Aop_\tau$ with $\tau=\frac{1}{n}$ for
    several values of $n$.  Solid lines are from sympy calculations
    and the dotted lines are from numerical power method in which the
    interval is divided into 3,000 pieces.}
  \label{fig:eigenfunctions}
\end{figure*}

\begin{figure*}
  \centering
  %\framebox[0.75\textwidth]{Plot here}
  \resizebox{0.75\textwidth}{!}{\includegraphics{eigenvalues.pdf}}
  \caption{A plot that illustrates dependence of $\lambda_{\tau}$ on
    $\tau$.}
  \label{fig:eigenvalues}
\end{figure*}

\subsection{Probabilities for Discrete Time}
\label{sec:discrete_t}

I begin by analyzing a process with a countable number of values for
the independent variable $t$ with uniform spacing $\tau$.  In this
analysis, I let $\lambda$ denote the largest eigenvalue.  By the Peron
Frobenius Theorem, the corresponding eigenfunctions can be chosen to
be non-negative with normalization that give the stationary
distribution for that process as
\begin{equation}
  \label{eq:stationary}
  \stationary(g) = \leftfunction(g) \rightfunction(g).
\end{equation}
The normalization condition is
\begin{equation*}
  \int \stationary(g) dg = 1.
\end{equation*}
The probability density at a particular sequence of $n+1$ samples from
$g_0$ to $g_n$ is
\begin{equation}
  \label{eq:joint}
  \density(g_0, g_1, \ldots , g_n) = \leftfunction(g_0) \prod_{k=1}^n
  \Aindicate(g_{k-1},g_k) \rightfunction(g_n).
\end{equation}
Thus the conditional density at $g_1$ given $g_0$ after $1$ step is
\begin{align}
  \density(g_1 | g_0) &= \frac{\leftfunction(g_0)
  \Aindicate(g_0,g_1) \rightfunction(g_1)}{ \int \leftfunction(g_0)
  \Aindicate(g_0,h) \rightfunction(h) dh} \nonumber \\
  \label{eq:one_step}
  &= \frac{\Aindicate(g_0,g_1)
    \rightfunction(g_1)}{\lambda \rightfunction(g_0)}.
\end{align}
Similarly, for two steps I find
\begin{align}
  \density(g_2 | g_0) &= \frac{\leftfunction(g_0) \left( \int
      \Aindicate(g_0,h_1) \Aindicate(h_1,g_2) dh_1 \right)
      \rightfunction(g_2) } 
      { \leftfunction(g_0) \int
      \Aindicate(g_0,h_1) \Aindicate(h_1,h_2) \rightfunction(h_2) dh_1
      dh_2} \nonumber \\
  &= \frac{ \left( \int
      \Aindicate(g_0,h_1) \Aindicate(h_1,g_2) dh_1 \right)
      \rightfunction(g_2) } 
      { \lambda^2  \rightfunction(g_0)} \nonumber \\
  \label{eq:two_step}
  &\equiv \frac{\delta(g_0) \Aop \Aindicate(\cdot,g_2)
    \rightfunction(g_2)}{\lambda^2 \rightfunction(g_0)},
\end{align}
and for $n$ steps
\begin{equation}
  \label{eq:n_step}
  \density(g_n | g_0) = \frac{\delta(g_0) \Aop^{n-1} \Aindicate(\cdot,g_n)
    \rightfunction(g_n)}{\lambda^n \rightfunction(g_0)}.  
\end{equation}
Note that
\begin{equation*}
  \delta(g_0) \Aop \Aindicate(\cdot,g_2) \equiv \int
  \Aindicate(g_0,h_1) \Aindicate(h_1,g_2) dh_1
\end{equation*}
is the length of the interval of $g_1$ values that are possible given
the values $g_0$ and $g_2$, and that the generalization
$\delta(g_0) \Aop^{n-1} \Aindicate(\cdot,g_n)$ is the volume in
$\REAL^{n-1}$ of tuples $(g_1,g_2,\ldots,g_{n-1})$ that are possible
given $g_0$ and $g_n$.

\subsection{Continuous Time}
\label{sec:continuous_t}

The principal task of this paper is to determine if a continuous time
limit of \eqref{eq:n_step} exits.  In particular, does the right hand
side of \eqref{eq:conditional} converge?
\begin{align}
  \label{eq:conditional}
  \density(g_t|g_0) = &\lim_{n \rightarrow \infty}  \frac{\delta(g_0)
    \left( \Aop_\tau \right)^{n-1} \Aindicate_\tau(\cdot,g_t)
    \rightfunction_\tau(g_t)} { \left( \lambda_\tau \right)^n
    \rightfunction_\tau(g_0)} \\
  & \text{ with } \tau=\frac{t}{n} \nonumber
\end{align}
Since Eqns.~\eqref{eq:limit_function} and \eqref{eq:limit_value}
provide\marginpar{Wrong}
$\lim_{\tau \rightarrow 0} \rightfunction_\tau(g) = e^{g-\frac{1}{2}}$
and $\lim_{\tau \rightarrow 0} \lambda_\tau = 1 - \frac{1}{e}$, all
that remains is to show that the limit in \eqref{eq:paths} exists.
\begin{align}
  \label{eq:paths}
  T_{t,g_0} & \equiv \lim_{n \rightarrow \infty} T_{t,g_0,n}, \text{
  where}\\
  T_{t,g_0,n} & \equiv \delta(g_0) \left( \frac{\Aop_{\frac{t}{n}}}
                {\lambda_{\frac{t}{n}}} \right)^n  \nonumber
\end{align}

For large $t$ the limit in \eqref{eq:paths} approaches a power method for
finding $\leftfunction$, and
\begin{equation*}
  \lim_{t \rightarrow \infty} T_{t,g_0} \propto \leftfunction.
\end{equation*}
However, if $t$ is smaller than $1 - g_0$ then the monotonic
constraint implies that
\begin{equation*}
  T_{t,g_0,n} = 0 ~\forall g,n: g > g_0 + t \text{ and } n > 0.
\end{equation*}
To simplify numerical calculations, I define
\begin{equation*}
  \tilde T_{t,g_0,n} \equiv \delta(g_0) \left( \frac{\Aop_{\frac{t}{n}}}
                {\lambda_0} \right)^n.
\end{equation*}
Numerical calculations of $\tilde T_{t=\frac{1}{4},g_0=\frac{1}{2},n}$
for several values of $n$ appear in Figure~\ref{fig:paths} and that
numerical evidence supports the conjecture that \eqref{eq:paths}
converges.
\begin{figure*}
  \centering
  \framebox[0.75\textwidth]{Plots here}
  %\resizebox{0.75\textwidth}{!}{\includegraphics{paths.pdf}}
  \caption{These plots of $\tilde T_{t=\frac{1}{4},g_0=\frac{1}{2},n}$
    for several values of $n$, suggest that $\lim_{n \rightarrow
      \infty}$ exists.}
  \label{fig:paths}
\end{figure*}


\newpage
\section{A simple example}
\label{sec:example}

Consider the adjacency matrix
\begin{equation}
  \label{eq:A}
  A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix},
\end{equation}
the schematic of the corresponding finite state Markov process in
Figure~\ref{fig:mt2}, and its probability matrix
\begin{equation}
  \label{eq:T}
  \T = \begin{bmatrix} a & b \\ c & 0 \end{bmatrix}.
\end{equation}

\begin{figure*}
  \centering
  \resizebox{0.4\textwidth}{!}{\input{mt2.pdf_t} }
  \caption{A schematic of the finite state Markov process specified by
  Equation~\eqref{eq:T}.}
  \label{fig:mt2}
\end{figure*}

Given values for $a$, $b$, and $c$, one can solve
\begin{equation}
  \label{eq:stationary_2}
  \mu = \mu \T
\end{equation}
for the stationary distribution $\mu$ and from there find the entropy
rate $h$ by
\begin{equation}
  \label{eq:rate}
  h = \sum_i \mu_i H(J|i)
\end{equation}
where
\begin{equation*}
  H(J|i) \equiv -\sum_j \T_{i,j} \log (\T_{i,j})
\end{equation*}
is the familiar $P\log(P)$ formula.

\subsection{Parameters of $\T$ that maximize $h$}
\label{sec:max}

For the maximum entropy values, given a trajectory length $N$, all
trajectories must have almost the same probability.  Since all
trajectories are composed of loops in a cycle basis, the following
equation must hold
\begin{equation}
  \label{eq:excycle}
  a^2 = b\cdot c.
\end{equation}
That combined with the normalization constraints for the two nodes
\begin{align*}
  a + b &= 1 \\
  c &= 1
\end{align*}
is sufficient to specify the solution
\begin{align*}
  a &= \frac{2}{1+\sqrt{5}} \approx          0.61803\\
  b &= \frac{\sqrt{5}-1}{1+\sqrt{5}} \approx 0.38197\\
  c &= 1.
\end{align*}
Solving \eqref{eq:stationary_2} for that specification of $\T$ yields
\begin{align*}
  \mu &=  \begin{bmatrix} \frac{1+\sqrt{5}}{2\sqrt{5}}, &
    \frac{\sqrt{5}-1}{2\sqrt{5}} \end{bmatrix} \\
  h &= \log\left(\frac{1+\sqrt{5}}{2}\right).
\end{align*}


\subsection{Topological entropy $\htop$}
\label{sec:htop}

The topological entropy of a directed graph specified by an adjacency
matrix $A$ is the rate at which the number of trajectories grows with
length.  If $n(t)\equiv\begin{bmatrix} n_1(t),&n_2(t)\end{bmatrix}$
denotes the number of trajectories of length $t$ that end in the two
states of Figure \ref{fig:mt2}, then
\begin{align*}
  n(1) &= \begin{bmatrix} 1,&1\end{bmatrix} \text{ and} \\
  n(t+1) &= n(t) A.
\end{align*}
Thus
\begin{align*}
  n(t) &= \begin{bmatrix} 1,&1\end{bmatrix} A^t \\
  \text{and }\lim_{t\rightarrow \infty} \frac{1}{t} \log n_1(t) &=
  \lambda
\end{align*}
where $\lambda$ is the largest eigenvalue of $A$.

Solving the eigenvalue problem yields
\begin{equation*}
  \htop = \log\left(\frac{1+\sqrt{5}}{2}\right).
\end{equation*}

The topological entropy provides one more equation of constraint on
the equations that the values in $\T$ must satisfy for maximizing $h$.
Rather than the single equation \eqref{eq:excycle}, the pair
\begin{equation*}
  - \htop = 2\log(a) = \log(b) + \log(c)
\end{equation*}
must hold.

\section{A Power Iteration Approach}
\label{sec:algorithm}

I will use the solution to this problem to characterize the
distribution of long trajectories in a system defined by an adjacency
matrix $A$.  I want the distribution to be uniform over all allowed
trajectories.

Given a list of allowed trajectories of length $2T$, I could
approximate the value of $\T_{i,j}$ by
\begin{equation*}
  \hat \T_{i,j} = \frac{_{2T}n_{T,T+1}(i,j)}{_{2T}n_{T}(i)},
\end{equation*}
where $_{2T}n_{T}(i)$ denotes the number of trajectories of length
$2T$ for which $i$ is the value at position $T$ and similarly
$_{2T}n_{T,T+1}(i,j)$ denotes the number of trajectories that have
values $i$ and $j$ at positions $T$ and $T+1$ respectively.  While the
exponential growth of the number of trajectories suggests that
implementing the estimate for large $T$ is not feasible, power
iterations for left and right eigenvectors of $A$ make it easy.

The number of allowed sequences of length say six that exactly matches
a given sequence is either zero or one and can be written as
\begin{equation*}
  _6n_{1,2,3,4,5,6}(a,b,c,d,e,f) = A_{a,b}A_{b,c}A_{c,d}A_{d,e}A_{e,f}.
\end{equation*}
The number of sequences that are only required to match at positions
three and four is
\begin{equation*}
 _6n_{3,4}(c,d) =  \sum_{a,b,e,f} {_6n}_{1,2,3,4,5,6}(a,b,c,d,e,f) =
 \sum_{a,f} \left(A^2\right)_{a,c} A_{c,d} \left(A^2\right)_{d,f}.
\end{equation*}
Similarly
\begin{equation*}
 _{202}n_{101,102}(c,d) = \sum_{a,f} \left(A^{100}\right)_{a,c} A_{c,d} \left(A^{100}\right)_{d,f}.
\end{equation*}

I define (and calculate) $R(t)$ and $L(t)$ recursively with the
following power iteration scheme
\begin{align*}
  L(1) &= \begin{bmatrix} 1,&1,&\cdots,&1,&1\end{bmatrix} \\
  N_L(t) &= \left| L(t) \right| \\
  L(t+1) &= \frac{L(t) A}{N_L(t)} \\
  R(1) &= L^{\text{T}}(1) \\
  N_R(t) &= \left| R(t) \right| \\
  R(t+1) &= \frac{A R(t)}{N_R(t)}.
\end{align*}
Note that as $t$ increases, $R(t)$ and $L(t)$ converge quickly to the
right and left eigenvectors respectively of $A$ that correspond to the
largest eigenvalue and I can calculate
\begin{equation}
  \label{eq:N}
   _{2(T+1)}n_{T+1,T+2}(c,d) \propto \left(L(T)\right)_c  A_{c,d}
   \left(R(T)\right)_d
\end{equation}
pretty easily.  With results from \eqref{eq:N} for a large enough $T$
to ensure convergence, one can calculate estimates $\hat \mu$ and
$\hat \T$ of the stationary distribution and transition probabilities
respectively as follows.
\begin{align*}
  \tilde L &= L(T) \\
  \tilde R &= R(T) \\
  \tilde P_{i,j} &= \tilde L_i A_{i,j} \tilde R_j \\
  \tilde m_i &= \sum_j \tilde P_{i,j} \\
  \hat \mu &= \frac{\tilde m}{\sum_i \tilde m_i} \\
  \hat \T_{i,j} &= \frac{\tilde P_{i,j}}{\tilde m_i}
\end{align*}
Maxentropic Markov chains appears in:\\
IEEE Transactions on Information Theory, \\
Date of Publication: Jul 1984\\
Author(s): Justesen, J.\\
Hoholdt, T.\\
Volume: 30 , Issue: 4\\
Page(s): 665 - 667

\section{Composition}
\label{sec:composition}

\begin{align*}
  \left( P(0\mapsto 1) \right)_{i,j} &= \frac{L_i A_{i,j} R_j}
    {\lambda z \frac{\sum_k L_i A_{i,k} R_k}{\lambda z}} \\
  &= \frac{A_{i,j} R_j}{\lambda R_i} \\
  \left( P(0\mapsto 2) \right)_{i,j} &=
     \sum_k \frac{A_{i,k} R_k}{\lambda R_i}
            \frac{A_{k,j} R_j}{\lambda R_k} \\
  &= \frac{(A^2)_{i,j} R_j}{\lambda^2 R_i}                                       
\end{align*}
If $B$ is the adjacency matrix for twice the step size of $A$,
sufficient conditions for consistency are
\begin{align*}
  B R &= \gamma R \\
  \frac{B_{i,j} R_j}{\gamma R_i} &= \frac{(A^2)_{i,j} R_j}{\lambda^2
                                   R_i} \text{or}\\
  (A^2)_{i,j} &= \frac{\lambda^2}{\gamma} B_{i,j}.
\end{align*}
The conditions are not plausible because if it is possible to get from
$i$ to $j$ in two steps, ie, $B_{i,j}=1$, the number of ways to get
from $i$ to $j$ in two steps will vary depending on the values of $i$
and $j$.

\section{Repetition}
\label{sec:repetition}

Here I've typed up hand written alternative notation.  For the one
step joint probability
\begin{align*}
  P_{j,k} &= \lim_{n\rightarrow \infty} \frac {\sum_{i,l} \left(A^n
            \right)_{i,j} A_{j,k}  \left(A^n \right)_{i,j} }{N_n}
            \text{ where the normalization } N_n \text{ makes} \\
  \sum_{j,k}P_{j,k} &= 1   \\
  P_{j,k} &= L_j A_{j,k} R_k \text{ where } L,R \text{ are the left and right
    eigenvectors of } A \\
  P_{k|j} &= \frac{L_j A_{j,k} R_k}{\sum_k L_j A_{j,k} R_k} \\
  &= \frac{L_j A_{j,k} R_k}{\lambda L_j R_j} \\
  &= \frac{A_{j,k} R_k}{\lambda R_j}
\end{align*}
For $n$ steps rather than one, the same calculations yield
\begin{align*}
  P_{k|j}(x(n)|x(0)) &= \frac{\left( A \right)^n_{j,k} R_k}{\lambda^n
                       R_j} \\
  &\propto \left( \delta_j \cdot A^n \right)_k R_k
\end{align*}


\end{document}

%%%---------------
%%% Local Variables:
%%% eval: (TeX-PDF-mode)
%%% End:
